load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

location_fc="/Users/ivana/WORK_BSC/data_to_publish/EC-Earth3P/control/"
location_fw="/Users/ivana/WORK_BSC/data_to_publish/EC-Earth3P/lowice/"

filtw = systemfunc ("ls "+location_fw+"toafluxes_allmons_ym_ensmean_y1130_tm.nc")
filtc = systemfunc ("ls "+location_fc+"toafluxes_allmons_ym_ensmean_y1130_tm.nc")

ftw    = addfiles (filtw, "r")
ftc    = addfiles (filtc, "r")

filsw = systemfunc ("ls "+location_fw+"surffluxes_allmons_ym_ensmean_y1130_tm.nc")
filsc = systemfunc ("ls "+location_fc+"surffluxes_allmons_ym_ensmean_y1130_tm.nc")

fsw    = addfiles (filsw, "r")
fsc    = addfiles (filsc, "r")

;print(fw)
;print(fc)
;*******************************************;
;Read in variables
;*******************************************
lat=fsc[0]->latitude
lon=fsc[0]->longitude
time=fsc[0] ->time
rad  = 4.0 * atan(1.0) / 180.

clat=dble2flt(cos(lat(:,0)*rad))

ntime= dimsizes(time)

rsdtc=ftc[:]->rsdt
rsdtw=ftw[:]->rsdt

rsutc=ftc[:]->rsut
rsutw=ftw[:]->rsut

rlutc=ftc[:]->rlut
rlutw=ftw[:]->rlut


rsdsc=fsc[:]->rsds
rsdsw=fsw[:]->rsds

rsusc=fsc[:]->rsus
rsusw=fsw[:]->rsus

rlusc=fsc[:]->rlus
rlusw=fsw[:]->rlus

rldsc=fsc[:]->rlds
rldsw=fsw[:]->rlds

hflsc=fsc[:]->hfls
hflsw=fsw[:]->hfls

hfssc=fsc[:]->hfss
hfssw=fsw[:]->hfss


toanetc=rsdtc(0,:,:)
toanetc=rsdtc(0,:,:)-rsutc(0,:,:)-rlutc(0,:,:)    ; positive down

surfnetc= rsusc(0,:,:)
surfnetc= rsusc(0,:,:)-rsdsc(0,:,:)+rlusc(0,:,:)-rldsc(0,:,:) +hfssc(0,:,:) +hflsc(0,:,:)  ;positive up

toanetw=rsdtw(0,:,:)
toanetw=rsdtw(0,:,:)-rsutw(0,:,:)-rlutw(0,:,:)   ; positive down

surfnetw= rsusw(0,:,:)
surfnetw= rsusw(0,:,:)-rsdsw(0,:,:)+rlusw(0,:,:)-rldsw(0,:,:) +hfssw(0,:,:) +hflsw(0,:,:)  ;positive up

tempsc=surfnetc
tempsc=surfnetc+toanetc

tempsw=surfnetw
tempsw=surfnetw+toanetw
printVarSummary(tempsc)
printVarSummary(tempsw)


;tempsc!0 = "lat"
;tempsc!1 = "lon"
;tempsw!0 = "lat"
;tempsw!1 = "lon"

tempsc_z=tempsc(:,0)
tempsc_z=dim_avg(tempsc(j|:,i|:))
tempsw_z=tempsw(:,0)
tempsw_z=dim_avg(tempsw(j|:,i|:))
;************************
; calculate the anomaly:
;************************

  finfo_ano=tempsw_z
  finfo_ano=tempsw_z - tempsc_z

  printVarSummary(finfo_ano)
  printVarSummary(clat)

  finfo_ano=finfo_ano*clat

  finfo_ano!0="lat"
  finfo_ano&lat=lat(:,0)
  ;print(finfo_ano)

;________________
;write data

;name ="columnnet_ano_zonal_ccsm4som"
;outf = addfile(name+".nc","c")
;outf->time = time
;outf->lat  = lat
;outf->finfo_ano   = finfo_ano

;****************************************
; Create plot
;****************************************
wks = gsn_open_wks("pdf","zonal_columnnetflux_ecearth3p")

   plot    = new (1,"graphic")

   res                      = True          ; individual plot
   res@gsnDraw              = False
   res@gsnFrame             = False
   res@xyLineThicknessF     = 2.0
   res@xyMonoLineColor         =True
   res@xyLineColor         = "black"
   res@gsnYRefLine           = 0.0             ; create a reference line
   res@gsnAboveYRefLineColor = "red"              ; above ref line fill red
   res@gsnBelowYRefLineColor = "blue"             ; below ref line fill blue

   res@xyDashPattern        = 0                  ; Make curves all solid

   res@vpHeightF= 0.4                    ; change aspect ratio of plot
   res@vpWidthF = 0.8

   res@trYMinF=-1
   res@trYMaxF=1.5

   res@tiXAxisFontHeightF = 0.010
   res@tiYAxisFontHeightF = 0.015

   res@gsnLeftString = ""
   res@tiYAxisString = "TOA net flux ano"
   res@tiXAxisString = ""

   plot(0)  = gsn_csm_xy (wks,lat(:,0),finfo_ano,res)

;***** make panel *****


   resP                     = True          ; panel resources
   resP@gsnMaximize         = True

   resP@txFontHeightF = 0.015
   gsn_panel(wks,plot,(/1,1/),resP)

